mcmc.plot <- function(file, act = "pairs", pthin = 1, useCols = NULL)
{

   var <- readList(paste(file,'.rep', sep=''))
   inData <- read.table(paste(file,'.mcmc.out', sep=''),col.names=var$mcnames)

    x = 1:nrow(inData)
    pthin = max(pthin, 1)
    xthin = seq(1, nrow(inData), pthin)
    if (pthin > 1)
        inData = inData[xthin, ]
    if (!is.null(useCols))
        inData = inData[, useCols]
    nc = ncol(inData)
    puce = "#cc8899"
    clrs = c("blue", "red", "green", "magenta", "navy", puce)
    clrs = rep(clrs, nc)[1:nc]
    resetGraph()
    panel = function(x, y) {
        len = length(x)
        points(x[2:len], y[2:len], pch = 21, col = "grey", bg = "gainsboro",
            cex = 0.8)
        points(x[1], y[1], pch = 21, bg = puce, cex = 1.2)
    }
    if (act == "pairs")
        pairs(inData, panel = panel, gap = 0)
    if (act == "eggs")
       plotFriedEggs(inData, eggs=TRUE, rings=TRUE, levs=c(0.01,0.1,0.5,0.75,0.95),
                   pepper=200, replace=FALSE, jitt=c(1,1), bw=25, histclr=NULL)
    if (act == "acf") {
        expandGraph(mfrow = c(1, 1), mar = c(3, 3.5, 0.5, 0.5))
        #clrs = rep(c("blue", "red"), nc)
        plotACF(inData, clrs = rep(clrs, nc)[1:nc],
            lwd = ifelse(nc > 6, 1, 2))
        mtext("Correlation", side = 2, line = 2.25, cex = 1.2)
        mtext("Lags", side = 1, line = 1.5, cex = 1)
    }
    if (act == "trace") {
        inData = cbind(x = xthin, inData)
        expandGraph(mfrow = c(nc, 1), mar = c(0, 0, 0, 0), oma = c(4,
            4.5, 0.5, 0.5))
        for (i in 1:nc) {
            plotTrace(inData[, c(1, i + 1)], clrs = clrs[i],
                xaxt = "n")
            axis(1, labels = par()$mfg[1] == par()$mfg[3])
            mtext(names(inData)[i + 1], side = 2, line = 3, cex = 1)
        }
        mtext("Valores secuenciales en la cadena de Markov", side = 1, outer = TRUE,
            line = 2, cex = 1)
    }
    if (act == "dens") {
        rc = PBSmodelling:::.findSquare(nc)
        expandGraph(mfrow = c(rc[1], rc[2]), mar = c(2, 2, 0,
            0), oma = c(1, 1.75, 0.5, 0.5), mgp = c(1.5, 0.2,
            0))
        for (i in 1:nc) {
            plotDens(inData[, i], clrs = clrs[i])
            addLabel(0.95, 0.95, names(inData)[i], adj = c(1,
                1), cex = 1.2)
        }
        mtext("Kernel Density", outer = TRUE, side = 2, line = 0.2,
            cex = 1.2)
        mtext("Parameter estimates", outer = TRUE, side = 1,
            line = -0.5, cex = 1)
    }
    invisible()
   return(inData)
}



read.fit<-function(file){
  # Function to read a basic fit
  ret<-list()
  parfile<-as.numeric(scan(paste(file,'.par', sep=''),
                      what='', n=16, quiet=TRUE)[c(6,11,16)])
  ret$nopar<-as.integer(parfile[1])
  ret$nlogl<-parfile[2]
  ret$maxgrad<-parfile[3]
  file<-paste(file,'.cor', sep='')
  lin<-readLines(file)
  ret$npar<-length(lin)-2
  ret$logDetHess<-as.numeric(strsplit(lin[1], '=')[[1]][2])
  sublin<-lapply(strsplit(lin[1:ret$npar+2], ' '),function(x)x[x!=''])
  ret$names<-unlist(lapply(sublin,function(x)x[2]))
  ret$est<-as.numeric(unlist(lapply(sublin,function(x)x[3])))
  ret$std<-as.numeric(unlist(lapply(sublin,function(x)x[4])))

  ret$cor<-matrix(NA, ret$npar, ret$npar)
  for(i in 1:ret$npar){
    ret$cor[1:i,i]<-as.numeric(unlist(lapply(sublin[i],
      function(x)x[5:(4+i)])))
    ret$cor[i,1:i]<-ret$cor[1:i,i]
  }
  ret$cov<-ret$cor*(ret$std%o%ret$std)
  return(ret)
}

matrix.ploteo <- function(mt, ...){
     x <- mt$cor
     min <- min(x)
     max <- max(x)
     yLabels <- as.character(mt$names)
     xLabels <- as.character(mt$names)
     title <-c()
  # check for additional function arguments
  if( length(list(...)) ){
    Lst <- list(...)
    if( !is.null(Lst$zlim) ){
       min <- Lst$zlim[1]
       max <- Lst$zlim[2]
    }
    if( !is.null(Lst$yLabels) ){
       yLabels <- c(Lst$yLabels)
    }
    if( !is.null(Lst$xLabels) ){
       xLabels <- c(Lst$xLabels)
    }
    if( !is.null(Lst$title) ){
       title <- Lst$title
    }
  }
# check for null values
if( is.null(xLabels) ){
   xLabels <- c(1:ncol(x))
}
if( is.null(yLabels) ){
   yLabels <- c(1:nrow(x))
}

layout(matrix(data=c(1,2), nrow=1, ncol=2), widths=c(4,1), heights=c(1,1))

 # Red and green range from 0 to 1 while Blue ranges from 1 to 0
 ColorRamp <- rgb( seq(0,1,length=256),  # Red
                   seq(0,1,length=256),  # Green
                   seq(1,0,length=256))  # Blue
 ColorLevels <- seq(min, max, length=length(ColorRamp))

 # Reverse Y axis
 reverse <- nrow(x) : 1
 yLabels <- yLabels[reverse]
 x <- x[reverse,]

 # Data Map
 par(mar = c(3,5,2.5,2))
 image(1:length(xLabels), 1:length(yLabels), t(x), col=ColorRamp, xlab="",
 ylab="", axes=FALSE, zlim=c(min,max))
 if( !is.null(title) ){
    title(main=title)
 }
axis(BELOW<-1, at=1:length(xLabels), labels=xLabels, cex.axis=0.7)
 axis(LEFT <-2, at=1:length(yLabels), labels=yLabels, las= HORIZONTAL<-1,
 cex.axis=0.7)

 # Color Scale
 par(mar = c(3,2.5,2.5,2))
 image(1, ColorLevels,
      matrix(data=ColorLevels, ncol=length(ColorLevels),nrow=1),
      col=ColorRamp,
      xlab="",ylab="",
      xaxt="n")

 layout(1)
}


read.myout <- function(file, start.year = NULL,...)
 {
     tt.chars <- scan(paste(file,'.rep', sep=''), sep = '', what='', quiet = TRUE)
     tt.nums <- as.numeric(tt.chars)
     list.length <- length(tt.nums[is.na(tt.nums)])
     pos <- which(is.na(tt.nums))
     object.counter <- 1
     out <- vector("list", list.length)

     while(object.counter <= list.length)
     {
         names(out)[object.counter] <- tt.chars[pos[object.counter]]
         if(object.counter < list.length)
         {
             out[object.counter] <- list(c(tt.nums[(pos[object.counter]+1):(pos[object.counter+1]-1)]))
         }   else
         {
             out[object.counter] <- list(c(tt.nums[(pos[object.counter]+1):length(tt.nums)]));
         }
         object.counter <-  object.counter + 1
     }

     if(is.null(start.year))
     {
         out$year <- c(seq(1,length(out[[1]])))
     } else
     {
         out$year <- c(seq(start.year,start.year+length(out[[2]])-1))
         out$yearpr <- c(seq(start.year,start.year+length(out[[2]])+out$years_projection-1))
     }
     out
 }


read.outfile<-function(path.string)
{
	tt.chars <- scan(path.string, sep = "", what = "")
	tt.nums <- as.numeric(tt.chars)
	list.length <- length(tt.nums[is.na(tt.nums)])
	out <- vector("list", list.length)
	object.counter <- 1
	element.counter <- 1
	while(object.counter <= list.length) {
		names(out[object.counter]) <- tt.chars[element.counter]
		print(tt.chars[element.counter])
		element.counter <- element.counter + 1
		tt.rows <- tt.nums[element.counter]
		element.counter <- element.counter + 1
		tt.cols <- tt.nums[element.counter]
		element.counter <- element.counter + 1
		if(tt.rows > 1) {
			out[[object.counter]] <- matrix(data = tt.nums[element.counter:(element.counter + tt.rows * tt.cols - 1)], nrow =
				tt.rows, ncol = tt.cols, byrow = T)
			element.counter <- element.counter + (tt.rows * tt.cols)
		}
		if(tt.rows == 1) {
			out[[object.counter]] <- c(tt.nums[element.counter:(element.counter + tt.cols - 1)])
			element.counter <- element.counter + tt.cols
		}
		object.counter <- object.counter + 1
	}
	out
}
